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Abstract 
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Obtainable computational efficiency is evaluated when using an Adaptive Mesh 
Refinement (AMR) strategy in time accurate simulations governed by sets of con- 
servation laws. For a variety of ID, 2D, and 3D hydro- and magnetohydrodynamic 
simulations, AMR is used in combination with several shock-capturing, conserva- 
tive discretization schemes. Solution accuracy and execution times are compared 
» ■ with static grid simulations at the corresponding high resolution and time spent on 

AMR overhead is reported. Our examples reach corresponding efficiencies of 5 to 
20 in multidimensional calculations and only 1.5 - 8 % overhead is observed. For 
AMR calculations of multi-dimensional magnetohydrodynamic problems, several 
strategies for controlling the V • B = constraint are examined. Three source term 
approaches suitable for cell-centered B representations are shown to be effective. 
For 2D and 3D calculations where a transition to a more globally turbulent state 
takes place, it is advocated to use an approximate Riemann solver based discretiza- 
tion at the highest allowed level(s), in combination with the robust Total Variation 
Diminishing Lax-Friedrichs method on the coarser levels. This level-dependent use 
of the spatial discretization acts as a computationally efficient, hybrid scheme. 
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1 Introduction 



Adaptive Mesh Refinement (AMR), an h- refinement method for solution- 
adaptive computations, is used in an ever increasing number of plasma physical 
and astrophysical simulations. The AMR algorithm and, specifically, all com- 
plications associated with ensuring global conservation properties when deal- 
ing with hyperbolic conservative systems, have been addressed in the work 
by Berger and Colella (0). In that pioneering paper, the authors stress that 
". . . it seems possible to implement a general code where the number of di- 
mensions is input". We make use of such a general code in all simulations 
presented here, where both the dimensionality and the particular system of 
conservation laws is selected in a pre-processing stage. The resulting versa- 
tility in the applications allows us to assess quantitatively how much can 
be gained from employing AMR in ID, 2D, or 3D hydro- and magnetohy- 
drodynamic (MHD) simulations, for a variety of spatial discretizations. Our 
coding effort benefits from the expertise gained through many general purpose 
and problem-tailored AMR software packages, which have been in active use 
and continued development for more than a decade. E.g., AMRCLAW ([29I ) by 
LeVeque and coworkers combines AMR with a selection of time-dependent hy- 
perbolic systems of partial differential equations, including both conservative 
and nonconservative systems. Berger and LeVeque (0) recently demonstrated 
the use of multi-dimensional wave propagation algorithms (|3ol: l28f) in 2D adap- 
tive computations on Cartesian and logically recta ngul ar, curvilinear grids. A 
related development effort by Walder and Folini (|49t l50r ) is AMRCART, a 
software package which intends to encompass 3D adaptive MHD scenarios. 
Dimension-independent library approaches for AMR simulations include the 
publicly available PARAMESH toolkit by MacNeice et al. 0), as well as C++ 
class framework approaches such as the AmrLib extension to the BoxLib li- 
brary |lo[). The code used here is the most recent extension to the Versatile 
Advection Code (see http://www.phys. uu.nl/ ^toth), briefly discussed in 



section 2.1. 

Our AMR implementation differs somewhat from the original algorithm in (0) 
and these differences will be mentioned in section 2.2. For several multidimen- 
sional calculations, we report on computational efficiency and solution accu- 
racy. We beneficially combine a fully upwind scheme on the highest AMR 
level with a Total Variation Diminishing Lax-Friedrichs discretization on all 
other levels. For 2D hydrodynamic and magnetohydrodynamic problems, we 
thereby typically gain a factor of 10 in execution time as compared to the 
equivalent high resolution static grid result. Our example simulations cover 
ID Riemann problems (Section 3.1), a 2D advection problem (Section 3.2.1), 
and both shock-dominated and more turbulently evolving hydro- and mag- 
netohydrodynamic multidimensional cases. Tables 3-4 collect our findings as 
presented in Section 3. 
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As our ultimate aim is to perform grid-adaptive, realistic 3D MHD simula- 
tions to investigate general plasma physical and astrophysical scenarios, we 
pay particular attention to strategies for controlling the solenoidal character 
of the magnetic field. In section 2.3, related efforts in AMR MHD simulations 
are discussed in relation to their specific treatment of this constraint. Recently, 
Balsara (2) and independently Toth and Roe (47) introduced a divergence-free 
reconstruction and prolongation strategy for multi-dimensional AMR MHD 
simulations which involves staggered magnetic field components. In this pa- 
per, we evaluate three alternative, simple approaches with cell-centered B 
components where V • B = is maintained to truncation error by the addi- 
tion of suitable source terms. Quantitative comparisons between AMR MHD 
calculations and high resolution static grid runs in section 3.2 demonstrate 
their effectiveness. 



2 Algorithm 



2.1 Versatile Advection Code, equations and solvers 



The Versatile Advection Code (VAC), initiated by Toth (|44j), is specifically 
designed for simulating dynamics governed by a system of (near-) conserva- 
tion laws. It's versatility resides in the choices offered in the geometry of the 
grid, the dimensionality, the physical application (including Euler and MHD 
equation modules), the spatial (|46| ) and temporal discretization (45[ 24), and 
the computer platform ([231) . 

VAC uses various second order shock-capturing numerical algorithms: the Flux 
Corrected Transport (FCT) scheme ([ill), the Lax-Wendroff type Total Varia- 
tion Diminishing (TVD) (Q) and the TVD-MUSClJII) schemes (TVDMU) 
with Roe-type approximate Riemann solvers (0; l39j). and the TVD Lax- 
Friedrichs (TVDLF) method (|52T). For the TVD type schemes different slope 
limiters are available, including the most robust minmod, and the sharper 
monotonized central (MC, also referred to as Woodward) limiter. For exact 
specifications of the algorithms, see (j46h and references therein. 

In multi-dimensional MHD simulations one needs to control the numerical 
value of V • B . VAC provides several algorithms, such as the eight-wave 
scheme (j36h . the projection scheme (0), and several types of constrained 
transport and central difference type discretizations, which were recently de- 
scribed and evaluated for static grid simulations by Toth (fiil ). In section 2.3, 
we discuss three strategies which easily carry over to an adaptive grid. Until 
now, VAC had no automated means for dynamic grid adaptation. Our new 
AMRVAC software package combines the versatility of VAC with the advan- 
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tages of adaptive mesh refinement. In this first implementation, only Cartesian 
grids and explicit time stepping is allowed. 

We demonstrate the flexibility of AMRVAC by using three equation modules 
in 1, 2 and 3 dimensions with different conservative high resolution spatial dis- 
cretization schemes. In the example simulations from Section 3, the equation 
modules are subsets of the system of magnetohydrodynamics, given by 



^ + V-(vp) = 5 p (1) 

do~v 

-|- + V • (vpv - B B ) + Vptota! = (2) 

dc 

_ + V-(ve + v Ptotal -BB -v) = 5 e (3) 

f)T2 

+ V .(vB -Bv) = S B . (4) 



These are expressed in the conservative variables density p, momentum density 
pv , total energy density e, and magnetic field B , with all additional sources 
and sinks collected in the source terms S p , S pv , S e and S#. In these equations, 
the total kinetic and magnetic pressure is written as 

Ptotal = (7 - 1) (e - ipv 2 - I B 2 ) + ^B 2 , (5) 

for an ideal gas with adiabatic index 7. We stress that for all simulations 
done, we time advance the relevant subset in the selected dimensionality only, 
which is chosen prior to compilation of the code. This is possible because the 
entire AMR algorithm is implemented in the LASY-syntax (|43h. making it 
effectively dimension- independent. Similar to VAC, the AMRVAC package is 
configured to dimensionality using the VAC preprocessor (a Perl script), prior 
to compilation. 

The time step is calculated to comply with the Courant-Friedrichs-Levy (CFL) 
stability criterion. If not stated otherwise, the Courant number is 0.9. The con- 
servative high resolution spatial discretization schemes used are FCT, TVDLF, 
TVD, and TVDMU. These shock-capturing schemes differ in the numerical 
representation of the fluxes Fj in the set of equations written as 

d t U + ]T <9 4 F,(U) = S(U, <9,U, didjV, x, t) . (6) 

i 

Here, the conservative variables are collected in U(x, t), and % and j run over 
1, 2 or 3 components of the spatial coordinate x. If not stated otherwise, 
multidimensional scenarios add fluxes from the different orthogonal directions 
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simultaneously, in contrast to a dimensionally split ([4 if ) approach. The source 
terms S can also be unsplit, i.e. added at the same time as the fluxes, or split, 
i.e. added before and after the update with fluxes. 



2.2 Adaptive Mesh Refinement 



Detailed descriptions of the AMR al gori thm can be found in Berger and 
Colella (0), Bell et al. (0), Friedel et al. (17), Steiner et al. (40), and many oth- 
ers. We summarize the essential steps and its user controlled features. Thereby, 
we pay specific attention to details which differ from the original Berger and 
Colella (|7|) approach. 



2.2.1 Time stepping and conservation 

In essence, AMR represents an automated procedure to generate or destruct 
- both controlled by the ensuing dynamics - hierarchically nested grids with 
subsequently finer mesh spacings Axi (index i — 1, . . . , D for a .D-dimensional 
problem). Up to a predefined maximal grid level Z max , consecutive levels I are 
characterized by even refinement ratios r/ with I 6 [2, Z max ], such that 

ri = Ax M-i = At *-i / 7 \ 
Ax ijt At/ 



To make the global conservation property computationally tractable, we stick 
to Cartesian grids and the 'proper nesting' criterion (0): (i) level I > 1 grid 
boundaries coincide with grid lines of / — 1 meshes; and (ii) we insist on the 
telescoping hierarchy of I — 2, I — 1 and level I grids, except at (non-periodic) 
physical boundaries. 

Consequently, two - related - issues for an AMR simulation are (1) how one 
time-advances on such a sequence of nested and thus overlapping grids, while 
(2) keeping the solutions consistent and conservative through a suitable 'up- 
date' process. The approach taken is illustrated in Figure 1, showing a hypo- 
thetical sequence of three time steps in a case where l max = 4. The schema 
is traversed from left to right, bottom to top, and with horizontal 'update' 
arrows preceeding vertical time 'advance' steps. Each time a level I > 1 has 
caught up in time with the underlying I — 1 grids, all level I — 1 cells that are 
covered by the finer grids are replaced by conservative averages. With nota- 
tion defined in Figure 2 for a 2D example where r\ = 2, this is achieved by 
simple averaging since e.g. the density in the underlying coarse cell must be 

p ij Ax li l- 1 Ax 2 ,l-l = J2mJ2n PmnAxi i lAx 2 ,l SO that fHj = Em En Pmn/ 'rf . Also, 

those / — 1 cells immediately adjacent to a level / grid boundary are modified 
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to comply with the interface fluxes used on the finer level. This step is crucial 
for restoring the global conservation property across the AMR grid hierarchy. 
Following the notation used in Figure 2, this involves a change in the temporal 
update for the density in coarse cell indexed by ij on level / — 1 bordered on 
the right hand side by level I cells. While originally 

fn+l) _ (n) A * ( F (n+\) At / («+£) p (n+§)\ 



the time centered numerical flux across the right cell edge F^t^ 
replaced by 



needs to be 



+ F V " V + F (n+ '^ + F (n+ ^ 

rn~^n m—^n+1 rn—^n rn—^n+1 



This flux fix restores the conservation. For more complete discussions of these 
update and flux fixing steps, we refer to (0). 

In our AMR implementation, each level solves the same set of equations but 
can exploit a different discretization scheme per level. This in contrast with the 
more general Adaptive Mesh and Algorithm Refinement (AMAR) approach 
introduced by Garcia et al. flih. who coupled a continuum description with a 
particle method on the finest level of the AMR hierarchy. AMAR is possible 
because the only communication between different level solutions is through 
(i) the update process and (ii) the filling of 'ghost' cells as used for imposing 
boundary conditions. 

The time advancing is further complicated by the automated regridding op- 
eration, as indicated by the grey circles in Figure 1. This regrid action is 
controlled by a simple criterion: when at least k time steps are taken on a 
certain level, this level is evaluated for refinement (k = 2 in the figure). This 
affects all higher level grids which may either suddenly appear (in Figure 1 af- 
ter the first and halfway in the second time step), disappear (after the second 
time step), or simply get rearranged or be left unchanged (halfway in the third 
time step). Of course, all level 1 grids remain unaltered, while the maximal 
allowed level Z max is never evaluated for further refinement. 



2.2.2 Refinement criteria and regridding 

Conceptually, regridding consists of two steps: in a first step, cells on level 
/ which need to be refined are identified - 'flagged' - through physically and 
numerically controlled criteria. The second step ensures that a properly nested 
new grid level I + 1 is formed which encompasses most, preferably all, of the 
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flagged cells. As pointed out above, the entire grid tree above a certain fixed 
base level /base m &y change in that process. With the current finest level in 
existence being Z fine , one descends in levels from I = min(7 fine , Z max — 1) down 

to / = /base- 



In the first step - the flagging procedure - a great variety of refinement cri- 
teria can be adopted. Powell et at (1371 ) use local measures of compressibility, 
rotationality, and current density in their 3D MHD simulations and a region 
is refined when chosen treshold values are exceeded. Ziegler (f53h uses simi- 
larly motivated, though more abstract, 'gradient value based' criteria where 
the slopes of functionals of the fluid variables, like the sound speed c s , control 
refinement. In cases where an (approximate) Riemann solver based method 
is exploited, one may flag cells in a 'wave- affected region'. The extent of this 
region can be computed from the wave strengths and speeds as calculated 
in the Riemann solver. Such criterion was succesfully exploited for unsteady 
Euler flows in (|14|). 



Our refinement criterion to 'flag' cells is a variant of the error estimation pro- 
cedure from Berger and Colella (J7J) , which is generally applicable to any system 
of equations and integration method used. The idea is to rely on Richardson 
extrapolation, by comparing a solution for a level I obtained by coarsening 
and integrating with one which is integrated first, and subsequently coars- 
ened. Bell et at (0) improved on this basic idea by additionally exploiting the 
uncoarsened solution vector as well, to avoid missing features that get lost 
in the averaging. Moreover, they allowed for a certain amount of user-forced 
refinement or derefinement. In our implementation, the integration method 
in the error estimation can be selected as any one-step, dimensionally un- 
split method, in combination with source terms. The criterion to flag cells can 
also exploit uncoarsened, high-resolution data, since we keep two full solution 
vectors U™ -1 and U™ on each level I < Z max , separated in time by the corre- 
sponding Atf~. For a fair variety of test problems (including the examples 
described below and those in Nool and Keppens (J33h. Keppens et at 
we found it sufficient to employ a first order version of the TVDLF method 
in the error estimation procedure, while integrating with a level-dependent 
mixture of TVD, TVDLF, or TVDMU on the full AMR hierarchy. The most 
crucial matter is to select the right combination of physics variables and mu- 
tual weighting factors to flag for. This is done as soon as a critical tolerance 
level etoi is exceeded. In the examples from section 3, we always list the chosen 
variables and tolerances. 



Given any procedure to identify cells which need refining, the AMR scheme 
must ultimately arrange them into properly nested 'rectangular' grids. Note 
that 'rectangle' is to be interpreted according to dimensionality. The basic 
steps to create this grid hierarchy were laid out in Berger and Colella (0) • In the 
original algorithm, each flagged cell was surrounded with a buffer zone of width 
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rebuff cells in each direction. The collection of flagged cells was subsequently 
made compatible with the nesting criterion. Some of the subtleties involved 
with ensuring this proper nesting are illustrated in Figure 3. Eventually, the 
collection of properly nested level I cells from which to form the level / + 1 
grids need to be processed into 'rectangles'. This latter step was based on a 
succession of bisections of the encompassing 'rectangle', followed by a merge 
step. Bisections were repeated till all grids had their ratio of flagged to total 
points within the rectangle greater than a predetermined efficiency E e [0, 1]. 
The merging was meant to balance the possible creation of many small grids. 
In its original form, this clustering algorithm created acceptable, but non- 
optimal new candidate grids which would often overlap and could end up 
with efficiencies below E, due to the merging process. Berger and Rigoutsos (J9|) 
designed a more sophisticated algorithm to overcome these problems, where 
the grid partitioning was done along edges where transitions from flagged to 
unflagged regions occured. The edges were detected using pattern recognition 
techniques exploiting signature arrays. This clustering algorithm has been used 
on realistic 3D problems by Bell et al. (0), Balsara ©, and various others. 
We chose to modify the original merge step from (0) to avoid the creation of 
overlapping grids at the same nesting level. At the same time, we introduced 
an enforced minimal efficiency E min 6 [0, 1] on all resulting new grids. It 
could be of interest to study how our modifications of the clustering algorithm 
compare with the Berger and Rigoutsos (0) approach on a variety of 2D and 
3D synthetic input clusters. In all simulations reported in what follows, we 
took E min = 0.5. 

The resulting new grid structure is then initialized from the solutions available 
on levels /base and up. For those cells which can not be copied directly from 
the same level on the original grid hierarchy, we use a conservative, minmod 
limited linear interpolation from the coarser level underneath. Due to the 
nesting property, this is always possible. In a notation similar to the one used 
in Figure 2, this limited linear interpolation in a 2D case would read 



Pmn Pij T . ^iP T l^jP 



where Aj,p^' denotes the minmod limited slope in direction % involving left- 
sided, right-sided, and centered slopes determined from [pi-ij, Pij, Pi+ij 1. Note 
that this is a conservative operation and second order accurate in smoothly 
varying density regions. 
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2.3 MulU- dimensional AMR MHD: the V • B =0 constraint 



Multidimensional MHD computations that make use of AMR have only been 
emerging over the last few years. Such computations must handle the addi- 
tional complication of maintaining a solenoidal magnetic field. Two-dimensional 
adaptive calculations for incompressible MHD flows were presented by Friedel 
et al. fl7h- Those authors borrowed projection method techniques from incom- 
pressible Navier-Stokes simulations (HI). Details of these projection methods 
and how to handle them on adaptive grid hierarchies have been discussed in 
Almgren et al. ((J). 



Steiner et al. (|40f l performed compressible MHD calculations using Flux Cor- 
rected Transport ((ill ) in two space dimensions. They combined AMR with a 
staggered representation where the magnetic field components are defined at 
cell interfaces while all other conserved quantities are cell-centered. Balsara (0) 
has presented and applied a divergence-free AMR MHD scheme where the 
same staggering was used. The key elements in his scheme are a divergence 
free prolongation strategy to transfer coarse solution information to finer grids, 
and a divergence free restriction step to update coarse grids from overlyin 
finer levels. Similar to the flux fix-up step inherent in the Berger and Colella 
method to ensure the global conservation property, a similar electric field cor- 
rection step is needed at mesh level interfaces to keep the overall solutions 
consistent. 



Powell et al. (|37l ) presented full 3D ideal MHD scenarios where solution- 
adaptive refinement is combined with an eight-wave Roe-type approximate 
Riemann solver method. Representative calculations of the magnetized solar 
wind impacting on a planetary magnetosphere have been reported by those au- 
thors. They advocated the use of an eight-wave formulation, which maintains 
the V • B = constraint to truncation error. The approach modifies the ap- 
proximate Riemann solver to 'propagate' VB errors with the plasma velocity 
and adds non-conservative source terms to the MHD equations proportional 
to V • B . The latter was found to work well for non-Riemann solver based 
methods as well (|46h. and is trivial to carry over to an adaptive grid. Recently, 
Janhunen (jiol ) and Dellar |l5h formulated various arguments to only add the 
V ■ B related source term to the induction equation, which restores momen- 
tum and energy conservation. Yet another approach has been forwarded by 
Marder (jl^), which relies on diffusing the numerically generated divergence 
at its maximal rate. Dedner et al. (jlfl ) revisited this 'parabolic' approach and 
introduced various variants for hyperbolic divergence cleaning. In summary, 
three 'source term' approaches to handle the V • B constraint have respective 
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sources (S p , S pv , S e , S B ) given by 



\ 



-(V-B)B 
-(V-B)B • v 
V "(V-B)v J 









pow 



v-(v-b) V ;. 








j an 



B • C d Ax 2 V (V • B ; 
^ C d Ax 2 V (V • B ) j 



dif 



In this paper, we apply these three simple approaches described by (8) in our 
cell-centered magnetic field representation. We cross-validate them on several 
2D MHD problems using AMR. The coefficient C d = 0.2 in all cases reported. 
This diffusive approach could also be altered to keep the energy equation fully 
conservative. 



3 Results and efficiency evaluation 



In this section, we present simulations covering ID to 3D evolutions. Our 
findings are discussed for each problem in what follows. Refinement is based on 
the density, if not stated otherwise. Tables 3-4 collect some relevant statistics, 
including execution times. All ID and 2D problems report timings done on a 
single 270 MHz MIPS R12000 CPU of an SGI Origin 200. The 3D cases have 
been obtained on single or multiple processors of the SGI Origin 3800, with 
500 Mhz MIPS R14000 CPUs. The efficiency of an AMR calculation is defined 
as the ratio between its execution time and the time needed to run the same 
problem on a static grid at the resolution of the highest level allowed Z max - 
We also report on the percentage of computing time spent on AMR-specific 
overhead. This overhead includes the total time spent on the combination of 
regridding operations (including the error estimation procedure and the steps 
needed for creating and initializing properly nested new grid hierarchies) and 
update operations, i.e. the replacement of underlying coarse cell values and 
the flux fixing of coarse cells neighbouring a finer grid level. This overhead 
measure leaves out two factors which are arguably AMR related, namely the 
time spent on storing fluxes across grid boundaries - which is negligible - 
and the time spent on enforcing boundary conditions by means of ghost cells. 
On AMR grid hierarchies, the latter involves spatio-temporal interpolation 
procedures from underlying coarser levels for filling internal boundary regions 
surrounding level / > 1 grids. 
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3.1 ID simulations 



3.1.1 Hydrodynamic shock tube 



Our first example considers a Riemann problem for the ID Euler system from 
Harten ()19l). With a ratio of specific heats 7 = 1.4, the left and right states 



arc 



(p,fW x ,e) L = (0.445, 0.311, 8.928) for x < 8 
(p,po x ,e) R = (0.5, 0, 1.4275) for x > 8 



in the x G [0, 14] domain with open (zero gradient) boundaries on both sides. 
Using 140 cells at level 1 and allowing for 5 levels with refinement ratios 
ri = 2, we strive to efficiently achieve the accuracy of a static grid simulation 
exploiting 2240 cells. In Figure 4, we show two solutions at t = 2: the top panel 
shows the density as obtained with FCT. In the bottom panel, the TVDMU 
scheme with MC/Woodward limiter was used on all levels. The extent of the 
level I > 2 grids at this time is indicated. We set the tolerance e to i = 0.005 
and regrid parameters k = 4, E = 0.6. 



As soon as the resulting rarefaction wave, contact discontinuity and shock 
are well-separated, higher level grids are only activated around the contact 
and the shock discontinuities. Locally, the solutions compare favorably with 
static, high-resolution results exploiting the same method: an equal amount 
of grid points resolving the contact and shock are found. We overplotted the 
high resolution static grid result in the bottom panel of Figure 4. Note that 
the FCT simulation does well in representing the discontinuities at the cost of 
introducing spurious oscillations, while the TVDMU scheme yields oscillation- 
free results. This may well be specific for the ETBFCT variant used here and 
described in Toth and Odstrcil (j46l). Similar conclusions were reached there 
comparing both methods for static grid simulations and in Steiner et al. (Jiol ) 
where AMR calculations used FCT. 



The execution times reported in Table 3 are for simulations up to time t = 
2, using a fixed time step At 5 = 0.000625. They demonstrate an obtained 
efficiency of 2.43 when using FCT, which improves to 2.95 when using TVDLF 
and up to 3.76 for the TVDMU scheme. This increase in AMR efficiency is 
reflected in a decrease of the AMR overhead from 25.7 % down to 19.8 %. 
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3.1.2 Magnetohydrodynamic Riemann problem 



Our second te st p roblem is the 1.5D MHD Riemann problem introduced by 
Brio and Wu (|13l). With fixed B x = 0.75 and 7 = 2, the constant states are 



(p, pv x , pv y ,e, B y ) L = (1, 0, 0, 1.78125, 1) for x < 0.5 
(p, pv x , pv y , e, B y ) R = (0.125, 0, 0, 0.88125, -1) for x > 0.5 

in the x G [0, 1] domain with open boundaries. We use 40 base cells and com- 
pare two means to achieve locally the resolution of a 1280 cell static grid. Using 
the TVDMU scheme and minmod limiting on all levels, we compare in Figure 5 
the density at t — 0.1 from a calculation with Z max = 6 and r; = 2 versus one 
where Z max = 4 and [r 2 , r 3 , r 4 ] = [4, 4, 2]. The tolerance was set to e to i = 0.001 
and regrid parameters k = 4, E = 0.8. Both solutions accurately capture the 
slow compound shock, the contact discontinuity, and the slow shock. Locally, 
the obtained solution essentially coincides with a high resolution static case, 
again overplotted in the bottom panel. The timings in Table 3 indicate only 
little differences in efficiency between the AMR simulations exploiting 4 and 
those exploiting 6 levels. This is true for the TVDLF and the TVDMU scheme. 
For the latter discretization, the AMR overhead drops to about 15 %. 

3.2 2D simulations 
3.2.1 Advection tests 

Our first 2D examples consider pure advection problems. On a doubly-periodic 
domain [0, 1] x [0, 1], a prescribed velocity field v = [u^i^] = [1, 1] diagonally 
advects a two-dimensional density profile. We perform two tests, one for a 
discontinuity dominated profile, and one for a smoothly varying density pulse. 

The discontinuous profile is a two-dimensional density field representing the 
VAC-logo, see Figure 6. Within the logo, the density is set to 2, while it is 0.5 
outside. We use a 50 x 50 level 1 grid, and allow for three levels with r\ = 2. 
The AMR control parameters are taken k — 4, E = 0.8, and e to i = 0.01. 

Figure 6 shows the initial t — grid structure, along with two snapshots at 
times t = 0.6 and t = 2, time-advanced with TVDLF and a MC/Woodward 
slope limiter. To demonstrate the maintained accuracy, cuts along the re- 
direction at y = 0.5 from times t = 0, 1, 2 are plotted as well. At each integer 
time, one should regain the t = situation. The numerical diffusion imme- 
diately smears the discontinuities, but the resulting density field is in exact 
agreement with a static grid result using 200 x 200 cells. As shown in Table 3, 
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the AMR simulation is twice as fast as a static grid run, with only 10 % of 
the execution time spent on regridding operations. 

For the timings reported using the FCT method, we used a dimensionally 
split approach. As pointed out by Bell et at (0), directional sweeping can be 
less favorable to the AMR calculations as additional boundary work is needed 
(extra top and bottom layers in the x-sweep must be taken along to start 
the y-sweep). The stencil of the FCT discretization for directionally unsplit 
multidimensional simulations requires filling the corner ghost cells, which we 
avoid with the split approach. 

Since the exact solution for this advection problem is known, we illustrate 
quantitatively that the AMR method reaches the same level of accuracy as 
a corresponding high resolution simulation. Table 1 summarizes the relative 
numerical errors obtained for a variety of AMR runs of different effective 
resolution. In general, we define the relative error for a variable u with respect 
to a reference solution w ref on a N x M grid as 



Su = 





Ui 


j-U 


ref 1 


2^i=i 2-,j=i 







(9) 



In this equation, the m t j values are to be obtained by extrapolating the AMR 
solution to a uniform grid of size N x M corresponding to its (highest level) 
effective resolution. When more than one variable is advanced, the average S 
is taken over all primitive variables u. 

Table 1 demonstrates that (1) the AMR runs achieve the same (first order) 
convergence behaviour on this discontinuity dominated problem as the static 
grid results; (2) the AMR runs are essentially identical to their high resolution 
static grid counterparts, with minute improvements noticable when refining 
more often (lowering k) and/or lowering the tolerance e to i- Note that it is 
advisable to take k consistent with the buffer zone width nwf used in the 
regridding process, as pointed out by Berger and LeVeque (J8j). As could be 
expected, the AMR efficiency is better for more grid levels than the 3-level 
one reported in Table 3: an efficiency of 6 is reached for the 800 x 800 AMR 
run exploiting 5 refinement levels. 

To demonstrate that we obtain second order accuracy on smooth solutions, 
we similarly advected a Gaussian bell profile, initially set to 

p(x, y, t = 0) = 2 exp (-100 \(x - 0.5) 2 + (y - 0.5) 



within a radius 0.2 from the center of the square domain and constant else- 
where. Again, we compare the obtained numerical solution after two full ad- 
vection cycles at time t = 2 with the exact solution. Table 2 shows that we get 
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second order convergence behaviour to the exact solution, on static and cor- 
responding AMR runs. The AMR simulations refined every second timestep 
(k = 2), while the tolerance was fixed at e to i = 0.001. At this tolerance, the 
5 level AMR result achieving an 800 x 800 effective resolution did not trig- 
ger level 5 grids on the entire bell profile. Note that we did obtain a solution 
with the specified accuracy 0.001, which shows the accuracy and efficiency 
of our error estimation procedure. Hence, this explains the apparent loss in 
convergence rate for this AMR run. 

Taken together with the ID test cases from above, these AMR convergence 
studies of 2D advection problems demonstrate some general trends. First, the 
effective resolution of the AMR grid is setting the overall accuracy of the ob- 
tained solution, as desired. Note that it is controlled by the combination of 
(1) the base grid resolution; (2) the number of refinement levels Z max ; and (3) 
the refinement ratio's in between consecutive levels 77. AMR runs of similar 
effective resolution which differ in the exact combinations taken of these three 
controlling elements should behave identical (cfr. the test problem shown in 
section 3.1.2). These conclusions are valid provided that the refinement crite- 
rion succesfully tracks all flow features. The latter is influenced by the choice 
of e to i and k. Table 1 demonstrates that only minor differences ensue when 
varying these parameters within reasonable bounds. A final set of AMR pa- 
rameters, consisting of the efficiencies E, E min , and ribuff plays a role in deter- 
mining the (highly varying) number and the structure of the grids on a certain 
level. There seems to be little reason to alter their values dramatically from 
run to run. 



3.2.2 Double Mach Reflection of a Shock 



A well-known shock dominated hydrodynamical test is the double Mach re- 



flection of a shock on a wedge, introduced by Woodward and Colella (|5ll). On 
a domain [0, 4] x [0, 1], a planar Mach 10 shock with post- and pre-shock states 
given by 

{p,pv x ,pv y ,e) post = (8, 8 x 8.25 sin 60°, -8 x 8.25 cos 60°, 563.5) 
(p,pv x ,pv y ,e) pie = (1.4, 0, 0, 2.5). 



makes an angle of 60° with the wedge, which is represented as a reflective 
boundary located at x G [1/6,4] and y — 0. The adiabatic index is 7 = 1.4. A 
self-similar pattern develops as the shock reflects off the wedge. The boundary 
conditions are fixed to the post-shock state at the left edge, open at the right 
edge, a partly fixed, partly reflective bottom, and a time- and space-dependent 
boundary value prescription along the top. 
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In the AMR simulation, we use TVDLF on a maximum of 4 levels, with an 
80 x 20 base resolution and refinement ratios r\ = 2. Furthermore, we took 
k = 6, E = 0.8, and a tolerance of e to i = 0.01 in the regrid process. Figure 7 
shows the density structure at time t = 0.2. We had to use a Courant number 
of 0.4 for the TVDLF discretization, since the x — 1/6 bottom boundary 
transition was not coincident with a cell edge at this resolution. As a result, 
a TVDMU simulation developed a carbuncle-type error at the location of the 
leading Mach shock intersecting with the wall. A parametric study of the 
AMR efficiency in terms of both execution time and memory needs for a large 
variety of effective grid resolutions is reported in Nool and Keppens (Jiih . 

Our timings indicate a significant gain in execution times - by a factor of 11 - 
between the corresponding 640 x 160 static grid simulation and the adaptive 
one. This high efficiency reflects the fact that the AMR process nicely succeeds 
in tracing the localized fronts, so that the domain coverage by level I > 1 grids 
remains low. In particular, level 4 grids initially cover only 5 % of the entire 
computational domain, with a modest increase to 15 % at the end of the 
calculation at t — 0.2. When raising the base resolution to 160 x 40 and again 
allowing for 4 levels with otherwise identical parameter settings, we reach an 
efficiency of 19.8 with still only 8% computing time spent on regridding. 



3.2.3 Hydrodynamic Rayleigh- Taylor Instability 

While the previous example represents a typical usage of AMR in hydro- 
dynamical shock governed evolutions, we now turn to a simulation where a 
transition to a fairly global, turbulent state takes place. In such cases, it is 
not a priori evident that we can benefit from the use of dynamically triggered 
mesh refinement. 

We simulate a 2D Rayleigh- Taylor instability on a [0, 1] x [0, 1] domain, with 
gravity g = — e y pointing downwards. As in Keppens and Toth (|22l). we start 
from a dense fluid with pdense = 1 above the interface yi nt = 0.8 + 0.05 sin27rx, 
resting on top of a light fluid with pn g ht = 0.1. Both fluids are at rest at 
t = and 7 = 5/3. The pressure field is set from a centered differenced 
hydrostatic balance dp/dy = —p, ensuring that the pressure about y int equals 
unity. This configuration is inherently unstable, while linear stability analysis 
in incompressible hydrodynamics predicts the shortest wavelengths growing 
the fastest. The staircased representation of the interface y- m t initiates such 
short wavelength perturbations, while its sinusoidal variation and the periodic 
side boundaries select a preferential longer wavelength variation. Top and 
bottom are impenetrable to the fluids. 

In the fixed 100 x 100 resolution and the (small) numerical diffusion 

associated with the TVDMU scheme suppressed the development of much 
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of the small-scale structure. Here, we repeat that simulation using AMR to 
locally achieve a resolution of 800 x 800. As we anticipate a need for higher 
resolution spreading across the entire computational domain, we want to make 
the cost associated with all intermediate levels I G [2, Z max — 1] as low as 
possible, without sacrificing accuracy on the highest level Z max . This can be 
done by using a level-dependent discretization scheme, such as the cheap but 
robust TVDLF method on all levels I < Z max — 1, in combination with a 
Riemann-solver based method, here TVDMU, on the highest level Z max = 5 
only. The base grid is 50 x 50 and all refine ratios are ri = 2. Other AMR 
parameters are as in the Woodward and Colella simulation. The snapshot 
of the logarithm of the density at time t — 1.9 shown in Figure 8 clearly 
demonstrates that we reached a very high effective Reynolds number regime. 

In Figure 9 (left panel) we plot the domain coverage for each grid level in 
the Rayleigh- Taylor simulation as a function of time. This figure confirms the 
visual impression of fine structure spreading across much of the computational 
volume, such that at time t = 1.9, more than 50 % of the area is at the high- 
est resolution allowed. Also, the level / = 2 grids cover the full domain at this 
time, making the / = 1 grid calculation completely obsolete. Clearly, since the 
physical process cascades to smaller scales, the strategy of minimizing compu- 
tational cost on lower lying levels through the use of a cheaper discretization 
pays off. In Table 3, we give an estimate of the execution time needed for a 
static 800 x 800 calculation from performing 20 time steps at this resolution. 
As the CFL-limited time step is observed to decrease monotonically in the 
AMR run, this allows to deduce both a lower and upper bound, and we report 
a mean value in the table. The corresponding efficiency of the AMR calcula- 
tion is about 8.3. Note also the low value of 3.9 % of adaptive grid related 
overhead. 



3.2.4 Orszag-Tang Vortex System 

The compressible evolution of the Orszag-Tang vortex system (jjil : l35h is a 
well-studied model problem where a shock-governed transition to an MHD 
turbulent state takes place. On a doubly periodic domain of size [0, 2n] x [0, 2tt], 
velocity "convection cells" described by v = [— sin y, sin x] are out of phase 
with magnetic islands resulting from B = [— siny, sin2x]. At time t = 0, the 
density and pressure take the uniform values p = 25/9 and p = 5/3 with 
7 = 5/3. We evolve this system till time t = 3.14, using four levels with a 
50 x 50 resolution at level 1. With r; = 2, we aim to find the 'reference' high 



resolution 400 x 400 solution shown in Fig. 16 from (42). As in that paper, 
our Figure 10 shows the temperature distribution, and the fine details agree 
closely. Other parameters used in the AMR simulation are k = 6, E = 0.7, 
and e to i = 0.025. 
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To demonstrate the advantages of the mixed strategy where TVDMU is used 
on the higher levels only, we report in Table 3 timings for 3 AMR simula- 
tions with / max = 4: one where we use TVDMU on all levels, one where we 
use it on levels / = 3 and / = 4, mixed with TVDLF on / = 1,2, and the 
advocated choice where TVDMU is active on level I = 4 alone. The efficiency 
increases from 5.85, to 6.28, to 7.38, respectively. In the right panel of Fig- 
ure 9, the coverage of the domain by the four grid levels is shown for the most 
efficient case. In agreement with the initially smooth evolution of the density 
(uniform at t — 0), higher level grids are only activated from about t = 0.5 
onwards. This represents an immediate gain in the obtained efficiency. As in 
the Rayleigh- Taylor simulation from above, the level I = 2 grids cover up the 
entire domain after some time. At the end of the simulation, this is about to 
happen for level 3 as well: a restart with a higher base level resolution should 
then be performed. 

This 2D MHD problem is also an excellent test for comparing the three differ- 
ent strategies (referred to as 'Powell', 'Janhunen' or 'diffusive' and defined by 
their respective sources in equation 8) for controlling V • B in AMR simula- 
tions. Since the exact solution to this problem is not available analytically, our 
best reference is to use the same high resolution solution used in Toth (J42J) for 
quantitative comparisons. As described in that paper, using a dimensionally 
split one-step TVD scheme with MC limiting, different approaches for con- 
trolling V • B still differ slightly at resolutions of 400 x 400. Static grid MHD 
simulations exploiting a projection scheme were found to yield almost al- 
ways the most accurate solutions. Repeating the Orszag-Tang simulation with 
4 AMR levels all exploiting the dimensionally split TVD scheme (for a fair 
comparison), the obtained relative errors as defined above at time t = 3.14 
with respect to a reference projection scheme solution are 



S pow = 0.0546 5 jan = 0.0808 S dii = 0.0465. (10) 



These values should be contrasted to remaining differences among two high 
resolution static grid simulations: between a constrained transport solution 
and a projection scheme approach, a relative error of 5 = 0.021 could be 
measured. Again, by lowering the tolerance e to i to 0.01, improvements up to 
<5dif — 0.02768 can be obtained. These values for the relative errors are also 
consistent with the values reported in Table 7 of Toth (42). Given the fact 
that for a considerable amount of time (for times t < 1), the AMR simulations 
could suffice with lower effective resolutions, the obtained relative errors are 
truly optimal. 
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3.2.5 Rotor 



A second ideal MHD problem frequently used for comparing static grid MHD 
simulations is the rotor problem introduced by Balsara and Spicer (J3|). We 
perform the 'first' rotor test as mentioned in Toth §4$) where the setup is as 
follows. On the domain [0, 1] x [0, 1], thermal pressure at t — is p — 1, and 
a uniform magnetic field has B x = and B y = 0. We set 7 = 1.4, and 

1 /2 

create a dense disk within r = [(x — 0.5) 2 + (y — 0.5) 2 ] < ro = 0.1 with 
p = 10 and v x = —2(y — 0.5)/r , v y = 2{x — 0.5)/r . For r > ri = 0.115 the 
fluid is at rest with p = 1. In between r < r < r±, linear profiles for density 
and angular speed connect both states. We run the simulation till t = 0.15. We 
strive for a 400 x 400 resolution, and use AMR exploiting 4 levels with k = 6, 
E = 0.7 and e to i = 0.025. Flagging is done on density and both magnetic field 
components with weighting factors 0.8, 0.1 and 0.1, respectively. A snapshot 
of the thermal pressure is shown in Figure 10. 

Since all important features are captured on the highest grid level at all times, 
adaptive simulations with 'source term' treatments to keep the magnetic field 
solenoidal yield effectively similar solutions as static runs. We quantify this 
latter statement as follows. We perform several AMR runs with a base grid of 
50 x 50 and a varying number of grid levels Z max = 2, 3, 4. In each case, we use 
TVDLF on the lowest level(s), with TVDMU and MC limiting on the highest 
level. The latter scheme was used for the reference 400 2 solution of Figure 18 
from §4$) . This static grid solution again used a projection scheme for V • B . 
As in the previous test case, even static grid results yield mutual differences of 
order 5 = 0.013 for a flux constrained transport versus a projection algorithm 
at this high resolution. When we then calculate the relative errors for AMR 
runs with effective resolutions increasing from 100 2 to 400 2 , these turn out to 
be 



^max 


= 2 


^pow 


= 0.0760 


^jan 


= 0.0794 


<5dif 


= 0.0753. 




= 3 


^pow 


= 0.0358 


$jan 


= 0.0412 




= 0.0373. 




= 4 


^pow 


= 0.0157 


$jan 


= 0.0248 




= 0.0107. 



These experiments demonstrate that each source term approach achieves first 
order convergence behaviour, as expected for this discontinuity dominated 
problem. Furthermore, the errors at effective resolutions of 400 2 are compara- 
ble to remaining differences between static grid solutions. As a final note, the 
advocated mixed discretization approach does not significantly sacrifice accu- 
racy by the use of the more diffusive scheme on lower levels. Indeed, when we 
repeat the simulation with / max = 4 and exploit TVDMU on all levels, the 
relative error for a diffusive approach becomes <5dif — 0.0096. From this and 
the previous test case, it seems that the diffusive approach is preferable over 
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both other source term treatments. 
3.2.6 MHD Kelvm-Helmholtz Instability 

Since the Orszag-Tang and the rotor system represented ideal MHD evolu- 
tions, we include a resistive MHD simulation taken from Keppens et al. (25). 
The initial condition has a uniform density p = 1 and pressure p — 1, a sheared 
and perturbed velocity field given by 

v = [0.645 tanh20y, 0.01 sin(2vrx)exp(-25y 2 )] 

and a discontinuously varying magnetic field B = ±0.129e x .. The magnetic 
field is thus of uniform strength, but changes sign at y = in the domain 
[0, 1] x [—1, 1]. With a resistivity rj = 10 -5 , high resolution simulations on 
static 400 x 800 grids demonstrated in (j^) how the nonlinear evolution is 
characterized by Kelvin-Helmholtz roll-up triggering tearing behavior. As the 
Kelvin-Helmholtz instability warps and amplifies the initial current sheet at 
y = 0, several magnetic islands develop due to local reconnection events, 
which are clearly visible in the density pattern (see Fig. 9 in reference (jUJ)). 
Shown in Figure 11 is a snapshot of the density pattern at t = 4.4 as obtained 
with an AMR simulation using 4 levels. Since it is vital to capture the initial 
dynamics at the y = interface accurately, we base the refinement criterion 
on the horizontal field component B x so that all levels are activated at t — 0. 
Other than that, parameters are k = 6, E = 0.8, e to i = 0.01. Instead of the 
one-step TVD method used in (j25h . here we use the combined TVDLF (on 
levels / < 3) and TVDMU (on I = 4) approach. The timing reported in Table 3 
is for a simulation up to t = 5, at which time the transition to a turbulent state 
has occured. Note the very high efficiency of 13.9 and the marginal regridding 
overhead. We used the Powell source terms in these calculations. 



3.3 3D Simulations 
3.3.1 3D Advection 

We made an efficiency assessment of a 3D advection problem by diagonally 
advecting a sphere of radius 0.2 in the unit cube where the density was set to 
2, while it was 0.5 external to the sphere. With triple periodic boundaries, the 
velocity v = (1,1,1) brings the sphere back to the center at t — 1. Similar 
to the 2D advection problem evaluated in Table 1, accuracy and convergence 
properties are equivalent to corresponding high resolution simulations. As re- 
ported in Table 4, a case where AMR exploits 3 levels, with a 40 3 base grid and 
refine ratios 2 and 4, demonstrates an efficiency of 19.9. The AMR overhead 
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has dropped to a rather negligible 4.14 %. This problem was run on a single 
CPU of the SGI Origin 3800. 



3.3.2 3D MHD problems 

As a first 3D MHD simulation, we recovered the evolution of a magnetized 
Rayleigh- Taylor instability as in Keppens and Toth Q). 3D AMR MHD 
simulations recovering previous high resolution Kelvin- Helmholtz unstable jet 
evolutions are documented in Nool and Keppens (jiih. In those simulations, 
the diffusive source terms for controlling V • B were exploited. The problem 
setup for the Rayleigh- Taylor study straightforwardly augments the 2D HD 
simulation in section 3.2.3 with a uniform horizontal magnetic field B = Q.le x 
and a 3D displacement of the interface separating the heavy from the light fluid 
according to y- 1Tl t = 0.8 + 0.05 sin 2irx sin 2irz. The magnetic field suppresses 
the formation of much of the fine scale structure. Here, the magnetic field 
constraint is handled by the Powell source method. 

Using AMR settings with k = 6, E = 0.7, e to i = 0.05, the timings for two 
AMR simulations with / max = 2 and / max = 3 recovering a 80 x 160 x 80 
static grid case are reported in Table 4. In both simulations, TVDMU is used 
on the highest level only, and a roughly fivefold efficiency is reached. Since 
the base grid is only 20 x 40 x 20, the problem setup immediately triggers 
level 2 grids on 20 % and 3 grids on 10 % of the total domain. These values 
increase to 50 and 20 %, respectively, at the end of the calculation. Using 
the average coverage values of 15 % and 35 %, and taking account of the 
relative refinement ratios, we can thus expect that the CPU time is at best 
reduced to 15 + 35/16 + 100/256 = 17.57 %, giving an optimal efficiency 
of 5.7, to be contrasted with the fivefold efficiency actually obtained. Since 
the runtime is dominated by the two highest levels, the obtained efficiency is 
truly optimal. The time spent on the regridding process is fully negligible: less 
than 2 % overhead was observed for this problem. Since one would typically 
use a higher base grid resolution, much better efficiencies can be expected to 
hold in practical 3D MHD simulations. Indeed, Ziegler (J53h demonstrated an 
efficiency of 45 on a purely kinematic 3D magnetic 'reconnection' problem 
(where only the induction equation (4) was solved) starting from a base grid 
of 100 x 120 x 60 while allowing for 3 grid levels. 

To demonstrate improved efficiencies at higher base grid resolutions, we set 
up a 3D MHD implosion problem on the cube [—0.5, 0.5] 3 with base grid 
resolution 60 3 . We again only allow 3 grid levels, with 77 = 2, so that an 
effective resolution of 240 3 is realized. At time t = 0, the problem setup has 
a uniform magnetic field B = ^3/5e y pervading a static medium where the 
pressure is 0.6 external to the sphere x 2 + y 2 + z 2 = 0.2 2 and 0.06 inside this 
sphere. The density is p = 1 except for a small off-center sphere within the 
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low pressure zone: where (x + 0.1) 2 + y 2 + z 2 < 0.04 2 the density is increased 
tenfold. We run this simulation using triple periodic boundaries up to time 
t = 0.4, corresponding to 125 level 1 discrete time steps At\. The solution 
at this time is shown in Figure 12, showing several quantities in cross-cuts 
at x — 0, y — 0, and z — 0. We used TVDLF on all levels and the diffusive 
source term treatment for the solenoidal constraint. Other AMR parameters 
were k — 2, E = 0.8, e to i = 0.01. This simulation was run on 16 processors 
of the SGI Origin 3800 and exploited OpenMP parallelization as explained in 



reference (|27f). 



During the entire run, the level 3 grids covered on average 7.3 % of the entire 
domain, reaching a maximum of 10.2 %, while the level I = 2 covered an 
average of 17.8 %, and a maximal coverage of 26.6 %. This already implies 
that this AMR run can at best achieve efficiencies in the range 8 — 11. For 
a fair comparison, we measured the wall clock time needed for performing 
one timestep at a fixed 240 3 resolution, while again exploiting 16 processors 
and using AMRVAC as a domain decompositioner. This high resolution static 
grid run needed then 78.9 seconds per time step. The average wall clock time 
needed for setting one Atx in the AMR calculation is 37.8 seconds. Taking 
into account the factor 4 in refinement from level 1 to level 3, the obtained 
efficiency is 8.3, exactly in its optimal range and indeed improved from the 
fivefold efficiency in the low base grid Rayleigh- Taylor problem. It is also 
worthwhile mentioning that the memory requirements for the 3-level AMR 
versus the static grid run are reduced by a factor of 6. Similarly, the space 
needed for storing one snapshot (in binary format where the three coordinates 
plus 8 unknowns per cell are saved in double precision) is a factor of 9 smaller 
in the AMR calculation. This also improves the time needed for IO operations: 
saving one snapshot from the AMR run on average took about 8 seconds, while 
that took roughly 51 seconds in the high resolution run. Hence, the overall 
advantages of performing grid-adaptive 3D MHD simulations are found in 
optimally reduced CPU times, drastically reduced memory requirements, and 
significantly smaller data volumes for eventual visualization purposes. 



4 Overview and outlook 



A dimension-independent implementation of the AMR algorithm, in combina- 
tion with the modular structure of the Versatile Advection Code, has allowed 
us to assess the gain in computing time for a variety of ID, 2D and 3D hydro- 
dynamic and MHD problems. The overhead associated with the automated 
regridding process is fully negligible for realistic 2D and 3D problems. A ten- 
fold efficiency is typical for the 2D cases shown, while the 3D timings indicated 
that very high efficiencies are indeed achievable: all 3D test problems demon- 
strated that the optimal efficiency reachable was in fact obtained. This will be 
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exploited in future AMR studies of correspondingly high resolution 3D MHD 
simulations of different, interacting fundamental plasma instabilities and stel- 
lar wind studies. It will be computationally beneficial to use Riemann solver 
based methods like TVDMU on the highest level only, in combination with 
the robust, more diffusive, but sufficiently accurate TVDLF discretization on 
lower levels. 



The AMR algorithm introduces a fair amount of control parameters, which 
do affect run strategies and the obtainable efficiency. As a general rule, two 
important parameter sets can be distinguished: (1) those controlling the ef- 
fective resolution (level I = 1 grid resolution, maximal allowed levels Z max , and 
consecutive refinement ratios rj); and (2) those determining the refinement 
criterion (the choice of variable(s) to refine on, the number k of time steps 
taken before regridding, and the tolerance in the criterion e to i). The physics 
problem at hand usually suggests suitable variables to use in the latter. Re- 
fining every second timestep k = 2 is a safe strategy, while the range for e to i 
in all problems studied above was e to i G [0.001,0.05], where the largest val- 
ues were used in the higher dimensional problems. As expected, the effective 
resolution plays a decisive role in the eventual efficiency. In cases where HD 
or MHD instabilities are simulated, their length scales will put constraints on 
the base grid resolution to use: it must be sufficiently high to allow for mode 
development. The localized nature of the problem will then help to determine 
what is a sufficient number of AMR levels l max : this may require some trial 
and error while analyzing the first few time steps to find an optimal value. 
Generally, raising the number of levels will beneficially influence the efficiency 



For mult i- dimensional MHD simulations exploiting AMR, three 'source term' 
approaches listed in equation (8) for controlling V • B to truncation error 
were investigated. Quantitative comparisons with reference high resolution 
simulations showed that all three work effectively for cell-centered B repre- 
sentations. For the tests included here, the diffusive approach was found to 
yield the smallest relative errors. As previously demonstrated by Toth |42T ) 
for static grid MHD simulations, we conclude that enforcing the solenoidal 
character of the magnetic field to machine precision in a particularly favoured 
discretization is not an absolute necessity for obtaining accurate multi-D MHD 
simulations exploiting AMR. The simple source term treatments offer a viable 
alternative to staggered field representations with considerable complications 
associated with constrained transport implementations on AMR grid hierar- 
chies (0). Projection scheme approaches would work accurately and reliably for 
AMR MHD simulations as well, but can be expected to be more complicated 
algorithmically and computationally costly. 
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Table 1 

Convergence study for the 2D advection of the VAC logo. The resolution indicates 
the effective resolution for AMR runs, refining every k (sub-)time steps. 



resolution Z max k e to i Static £> ex act 

100 2 1 0.12866 

2 4 0.01 0.000179 0.12867 

2 2 0.01 0.000215 0.12867 

2 2 0.005 0.000132 0.12866 
200 2 1 0.07646 

3 4 0.01 0.000301 0.07656 
3 2 0.01 0.000172 0.07648 

3 2 0.005 0.000102 0.07648 
400 2 1 0.04612 

4 4 0.01 0.001018 0.04676 
4 2 0.01 0.000222 0.04618 

4 2 0.005 0.000129 0.04617 
800 2 1 0.02781 

5 2 0.01 0.000506 0.02811 
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Table 2 



Convergence study for the 2D advection of the Gaussian bell profile, 
indicates the effective resolution for AMR runs (Z max > 1), refinin 
(sub-)time steps, at a tolerance e to i = 0.001. 


The resolution 
g every k = 2 


resolution 


^max 


<5cxact 


100 2 


1 


0.0238102 




2 


0.0238104 


200 2 


1 


0.0058374 




3 


0.0058379 


400 2 


1 


0.001548 




4 


0.001571 


800 2 


1 


0.000422 




5 


0.001271 
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Table 3 



Overview of all ID and 2D Simulations and Timin; 


gs (in 


seconds) . 






problem 
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2D HD Shock reflection 
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1372.82 


6.28 


3.49 




TVDLF-MU (3-1) 


4 


1167.91 


7.38 


3.94 


2D MHD Kelvin-Helmholtz 


TVDMU 


1 


173572 






400 x 800 x [0, 5] 


TVDLF-MU (3-1) 


4 


12477 


13.91 


2.36 
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Table 4 



Overview of 3D Simulations and Timings (in seconds). Note that the 3D Implo- 
sion problem mentions wall clock times for OpenMP parallelized code running on 16 
processors of the SGI Origin 3800. 

problem method l max timing eff. 


3D Advection 


TVDLF 


1 


172197 




320 3 x [0, 1] 




3 


8653 


19.90 


3D MHD 


TVDMU 


1 


218650 




Rayleigh- Taylor 


TVDLF-MU (1-1) 


2 


46659 


4.69 


80 x 160 x 80 x [0, 1] 


TVDLF-MU (2-1) 


3 


40726 


5.37 


3D MHD Implosion 


TVDLF 


1 


39437 




240 3 x [0, 0.4] 




3 


4733 


8.3 



Fig. 1. A hypothetical sequence of three time steps in an AMR simulation allowing 
for 4 nested refinement levels. Vertical arrows indicate a temporal advance, horizon- 
tal arrows correspond to update operations ensuring conservation and consistency, 
while the automated regridding operation is called at the position of the grey circles. 
Regridding is invoked when at least k = 2 time steps are taken. 



Fig. 2. Illustration of the update and flux fix step. An underlying coarse (level I — 1) 
cell value is merely replaced by a conservative average. A coarse cell bordering a 
finer level I grid needs its numerical flux across the edge x i+ i replaced by the more 
accurate fluxes used on the level I. 

Fig. 3. The regridding operation illustrated in 2D. Top left: 5 cells are flagged for 
refinement on 2 adjacent level I grids. The edge of the computational domain is 
indicated by thick solid lines. Each flagged cell is surrounded by a 2-cell buffer zone. 
Top right: Buffer cells not located on level I grids are discarded. A previously formed, 
enlarged level I + 2 grid (small rectangle) is to be projected for proper nesting in 
the level I + 1 grids being formed. Bottom left: discard cells which would violate the 
proper nesting criterion. Note the flagged cell which is lost! Bottom right: A possible 
grid structure for level I + 1, whereby two extra level / cells (hashed squares) are 
taken along. Note that the level I + 1 grids may overlap multiple level I grids. 

Fig. 4. The density structure at time t = 2 for the ID hydrodynamic Harten shock 
tube. Using AMR with 5 grid levels, we confront a solution obtained with FCT (top 
panel) with a TVDMU result (bottom panel) . The extent of the level I > 1 grids is 
indicated by the solid bars. A high resolution static grid result is overplotted as a 
solid line in the bottom panel. 

Fig. 5. The density variation at t = 0.1 in the Brio-Wu shocktube problem. A 
solution using 6 refinement levels (top panel) is compared with one where 4 levels 
and correspondingly higher refinement ratios were used (bottom panel). The extent 
of the level I > 1 grids is indicated by the solid bars and the bottom panel shows a 
high resolution reference result. 
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Fig. 6. Diagonal advection of the VAC-logo on a doubly periodic domain. We show 
the initial 3-level grid structure at top left, and snapshots at t = 0.6 (bottom 
left) and after two full advection cycles (bottom right). Cuts at y = 0.5 for times 
t = 0, 1, 2 (top right) demonstrate the maintained accuracy. 



Fig. 7. The density structure at t = 0.2 for the Woodward and Colella reflected 
shock problem, obtained with an AMR simulation using 4 grid levels and a 80 x 20 
base grid. 



Fig. 8. The logarithm of the density at time t = 1.9 shows the fine scale mixing 
process of a heavy fluid (dark) into a lighter one underneath under the influence of 
gravity. Only the level I = 5 grids are indicated. 



Fig. 9. For the Rayleigh- Taylor simulation from Figure 8 (left) and for the 
Orszag-Tang vortex simulation from Figure 10 (right), the domain coverage of the 
grid levels used in the AMR simulation as a function of time. 



Fig. 10. Left: the temperature at t = 3.14 for the Orszag-Tang magnetized vortex 
system as obtained with a 4 level AMR simulation. Right: the thermal pressure at 
t = 0.15 for the rotor problem as obtained with a 4 level AMR simulation. 



Fig. 11. The density structure at time t = 4.4 for a Kelvin-Helmholtz unstable 
shear flow configuration containing a current sheet shows the triggering of magnetic 
islands through tearing instabilities. Only the central part of the computational 
domain is shown. 



Fig. 12. 3D MHD implosion problem at time t = 0.4. Top left: in the plane z = 0, 
we show a schlieren plot of the pressure, the magnetic field lines, and the velocity 
field. Top right: schlieren plot of the density in y = 0. Bottom panels: density 
(left) and thermal pressure (right) variation in the plane x = 0. 
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